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Disordered fibrous networks are ubiquitous in 
nature as major structural components of living 
cells and tissues. The mechanical stability of net¬ 
works generally depends on the degree of connec¬ 
tivity: only when the average number of connec¬ 
tions between nodes exceeds the isostatic thresh¬ 
old are networks stable (ll|. Upon increasing the 
connectivity through this point, such networks 
undergo a mechanical phase transition from a 
floppy to a rigid phase. However, even sub¬ 
isostatic networks become rigid when subjected 
to sufficiently large deformations. To study this 
strain-controlled transition, we perform a com¬ 
bination of computational modeling of fibre net¬ 
works and experiments on networks of type I col¬ 
lagen fibers, which are crucial for the integrity 
of biological tissues. We show theoretically that 
the development of rigidity is characterized by a 
strain-controlled continuous phase transition with 
signatures of criticality. Our experiments demon¬ 
strate mechanical properties consistent with our 
model, including the predicted critical exponents. 
We show that the nonlinear mechanics of colla¬ 
gen networks can be quantitatively captured by 
the predictions of scaling theory for the strain- 
controlled critical behavior over a wide range of 
network concentrations and strains up to failure 
of the material. 

As shown by Maxwell, networks with only central-force 
interactions exhibit a transition from a floppy to rigid 
phase at the isostatic point, where the local coordination 
number, or connectivity ( 2 :) equals the threshold value of 
(z) = 2d in d dimensions [l|. At this point, the number 
of degrees of freedom is just balanced by the number of 
constraints, and the system is marginally stable to small 
deformations. The jamming transition l2|-l6| in granular 
materials and rigidity percolation |7Hld| in disordered 
spring networks are examples of such a transition. An 
important feature of these systems is the order of the 
transition. Jamming exhibits signatures of both first- and 
second-order transitions, with discontinuous behaviour of 
the bulk modulus and continuous variation of the shear 
modulus Ullllli. For networks of springs or fibres, the 
transition from floppy to rigid is a continuouy)hase tran¬ 
sition, in both bulk and shear moduli jl, S, 0, 0, 0, [l3|- 

Interestingly, the structural networks in biology almost 
always have connectivities below the central-force iso- 



Figure 1 : At zero strain, 7 , networks undergo a continu¬ 
ous transition from floppy to rigid at the isostatic threshold 
2 = 2 c- This connectivity threshold shifts to lower values for 
networks subject to shear strain 7 . This threshold defines a 
line 7 c( 2 :) of continuous transitions. We study here strain- 
induced transitions indicated by the vertical dashed line for 
2 well below 2 c. The insets show SEM (Scanning Electron 
Microscope) images of reconstituted collagen networks indi¬ 
cating points of 3-fold and 4-fold connectivities. The scale 
bar corresponds to 200 nm. 

static point. Networks such as the extracellular matrices 
of collagen that make up tissue are strictly sub-isostatic 
with respect to central forces: their typical connectivity is 
between 3 (local branching) and 4 (binary crosslinking), 
placing them below both 2D and 3D isostatic thresh¬ 
olds |l7l. Il8|. Such sub-isostatic networks can, neverthe¬ 
less, become rigid as a result of other mechanical con¬ 
straints, such as fibre bending 0, 0, [lH, or when sub¬ 
jected to external strain 0. The threshold strain, at 
which the transition occurs, depends on the nature of 
applied deformation (shear or tensile) and on the aver¬ 
age connectivity of the network, in particular, and other 
properties of its structure (Ij. However, the order of this 
strain-induced transition remains unclear. 

Here, we study the transition from floppy to rigid 
states of disordered sub-isostatic networks under simple 
shear. We show that these networks exhibit a line of sec¬ 
ond order transitions (see Fig. [1]) at a strain threshold 
7 c( 2 :), for connectivities 2 well below the isostatic thresh¬ 
old. Moreover, we demonstrate critical behaviour along 
this line, specifically in the scaling properties of the me¬ 
chanics, as well finite-size effects that reflect the under¬ 
lying divergent correlation length. To test the relevance 
of these predictions for real materials, we perform exper- 
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Figure 2: (a) Stiffness K in units of jji/l vs. strain 7 for 2D phantom triangular lattices with {z) ^ 3.2. The red dashed curve, 
starting from 7 = 7 c, is a sketch of the stiffness of a sub-isostatic network with ^ = 0 while the red solid curve is the affine 
limit with k = 00 . The inset shows the ^-dependence of the linear modulus Go with the dashed line indicating a unit slope, 
(b) Collapse of stiffening curves for different k. Black curve: the data set from (a) collapsed according to Eq. ([!]),/ = 0.75 
and (j) = 2.1. Blue curve: 2D Mikado network with {z) ^ 3.6, / = 0.84 and 0 = 2.2. Red curve: 3D phantom FCC lattice with 
[z) ~ 3.2, / = 0.8 and 0 = 2.2. In 3D, K is measured in units of fi/l^. The solid line has unit slope and the dashed line has 
slope f /(j) — 0.38. The inset shows the data collapse according to Eq. m for a disordered honeycomb lattice with (z) ~ 2.4 
with critical exponents / = 0.48 and (j) = 2.2. 


iments on reconstituted networks of collagen, the most 
prevalent protein in mammals and the mechanical basis 
of most tissues (lH. Although collagen has been widely 
studied for many years, the mechanical properties of col¬ 
lagen matrices remain poorly understood. We find that 
collagen networks show evidence of critical behaviour in 
their mechanical response to strain. Strikingly, the mea¬ 
sured shear modulus of these networks is in quantitative 
agreement with the critical behaviour of our model, in¬ 
cluding the predicted non mean-field critical exponents. 

We study computational models (See Methods) of fibre 
networks, based on both 2D and 3D lattice-based struc¬ 
tures miiii and Mikado networks in 2 D Hi. All net¬ 
works are, by construction, sub-isostatic and floppy in the 
absence of bending interactions j3| • The filaments have a 
stretching modulus, /i, and bending modulus, k. These 
parameters define a dimensionless rigidity k = 
where I is the lattice spacing (mesh size) in lattice-based 
(Mikado) networks. The networks are subjected to sim¬ 
ple shear strain 7 and allowed to relax by minimization 
of the total elastic energy per unit volume, 77, which is 
calculated using a discrete form of the extensible worm¬ 
like chain Hamiltonian [^. The stress and stiffness, in 
units of are obtained from 77 by cr = ^ and 

TC = respectively. 

In Fig. 0 we show the network stiffness K vs. strain 7 
of a triangular network with {z) 3.2 for different values 

of k. As sketched in Fig.[TJ these networks are character¬ 
ized by a continuous transition at a strain threshold 7 c, 


which is indicated in Fig. [2^ by the vertical dashed red 
line, above which the stiffness K increases continuously 
from zero for k = 0. This curve is approached for sys¬ 
tems with finite but decreasing as can be seen by the 
lower sets of black data points. This second-order phase 
behaviour is qualitatively analogous to the onset of ferro¬ 
magnetism on lowering the temperature below the Curie 
temperature, where the addition of a magnetic field re¬ 
sults in finite magnetization in the paramagnetic phase. 
More precisely, an additional energy Q such as the elas¬ 
tic bending stiffness of fibres j 3 | with a finite coupling 
constant k can stabilize otherwise floppy networks. As 
we show below, the stabilizing effect of k can be used 
to reveal the critical behaviour for sub-isostatic systems 
with z < Zc dit strains 7 7 c(^)- 

In the absence of the stabilizing effect of bending, i.e., 
for ^ = 0 , the continuous nature of the transition in 7 is 
apparent (Supplementary Fig. S2) in the critical scaling 
of the network stiffness K ~ in the regime where 

A 7 = 7 — 7^ > 0. For 7 < 7c, the effect of stabiliza¬ 
tion by bending leads to K k. These regimes can be 
summarized by the scaling form 

i^oc|A7Ka±(7|A7l^), (1) 

where Q± is a scaling function with the positive and neg¬ 
ative branches corresponding to A 7 > 0 and A 7 < 0, 
respectively. This scaling is analogous to that for the 
conductivity of random resistor networks and fibre net¬ 
works as a function of connectivity 0 , 123 . In Fig. H), 
we test this scaling relation by plotting vs. 
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c|A 7 | ^ [mg/mL] 


Figure 3: (a) Nonlinear stiffness vs. strain measured for collagen networks prepared at different concentrations (see legend). 
The inflection point ( 7 c) is in each case marked in red enlarged symbols. The dashed lines are the prediction of Eq. ilgTali . 
The fit values of k based on Eq. (ISBli are shown in the inset. The linear scaling of k with the concentration c is consistent 
with the predictions of our model (Supplementary Information). Also shown in the inset is the linear modulus Gq scaling with 
concentration as c^‘^. (b) Collapse of experimental stiffening curves for different concentrations of collagen (see legend). Solid 
line has slope 1.2 and the dashed line //0 = 0.36. Inset shows the weak dependence of critical strain on collagen concentration. 
The dashed line represents the model prediction of 7c ~ 


^|A 7 |“'^, according to Eq. ([T]). For x <C 1, Q+{x) is 
approximately constant and Q-{x) oc x. Since K must 
be finite at A 7 = 0, we also expect K ^ ^^^^^ 
as one observes from the critical branch in Fig. [Jb, con¬ 
sistent with Eq. ©• To show the generality of this re¬ 
sult, we also show in Fig. O) the corresponding data ob¬ 
tained from Mikado networks, as well as FCC-based 3D 
lattices. Strikingly, the data collapse with similar expo¬ 
nents / = 0.8 ± 0.05 and 0 = 2.1 ± 0.1. The average 
connectivity for the three different networks is chosen to 
be in the range 3.2-3. 6 , comparable to typical biopoly¬ 
mer networks such as collagen 

To test these predictions for networks of collagen type 
I, we measure the stiffness vs strain over a range of 
concentrations using shear rheometry (See Methods) as 
shown in Fig. [3^1. We also measured the average coordi¬ 
nation number of these networks, z = 3.3T0.1, consistent 
both with previous studies 0 and the values used above 
in our model. The parameter k in our model is naturally 
related to the protein concentration c (See Fig. [3^ and 
Supplementary Information) such that by rescaling our 
experimental K by the concentration we can compare 
with the simulations. In both experiment and simulation, 
we obtained the critical strain 7 c as the inflection point 
of the \ogK vs. log 7 curves (Supplementary Fig. Sla). 
This strain also coincides with the strain at which the 
non-afiine fluctuations in the network diverge (Supple¬ 
mentary Fig. Sib). Based on these considerations, our 
model predicts that the experimental stiffness K should 
be governed by the scaling relation 

i^/coc|A7|/a±(c/|A7|q. (2) 


In Fig. [ 3 ) 3 , we test this prediction by plotting iFjAyl ^ jc 
vs. cjAyl”"^. We find an excellent collapse with our pre¬ 
dictions for 3D networks shown in Fig. [2)3, assuming ex¬ 
ponents / = 0.8 and 0 = 2.1. 

In the limit of 7 ^ 0, using Q-{x) ^ x (Eq. ([1])), 
we obtain a scaling relation for the linear modulus and 
the concentration: Kjc ^ Since 7 c is expected 

to depend on the average connectivity, which may vary 
with concentration, we expect 7 c to show a possible 
concentration dependence. In fact, as seen in the in¬ 
set of Fig. [3)3, the experimental 7 c does show a weak 
7 c ^ dependence on the concentration. More¬ 

over, this dependence is consistent with the observed 
difference between the experimental and theoretical Q- 
branches of Fig. and Kjc^ is consistent with 
Kjc ^ and the exponents / and (p in Fig. [3)3, 

since jc ^ c^-^/D- 0 ) ^ 

To quantitatively capture the mechanics of fibrous net¬ 
works over the entire strain-range, one needs the scaling 
function Q±(x). As for ferromagnetic systems, one can 
obtain an approximate scaling function by numerical in¬ 
version of the equation of state @ (Supplementary Infor¬ 
mation), which in our case is given by 

k K f 

Here, ± corresponds to the two branches of the scaling 
function. In Fig. [3b, we show the predicted K vs. 7 ac¬ 
cording to Eq. ( 1 ^ . with k as the only fitting parameter 
(Supplementary Information). 


( 0 -/) 


(3) 
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Figure 4: Demonstration of the continuous transition of elas¬ 
ticity in a sub-isostatic network. Collapse of stiffening curves 
according to Eq. (p7|l for different system sizes W. Inset shows 
the scaling of | 7 c(oo) — 7 c(kF)| with system size. The dashed 
line has a slope of — 1 / 2 . 

We further test the critical behaviour by performing 
finite-size scaling, which is sensitive to the divergence of 
the correlation length. The modulus should follow the 
scaling relation (Supplementary Information) 

K oc W-^l''T± (|, (4) 

where W is the system size and J^± is a scaling function 
with the positive and negative branches corresponding 
to A 7 > 0 and A 7 < 0, respectively. In Fig. [S7l we 
plot iFFF^/^ vs. jAylVF^/^, showing consistency with 
Eq. ([S7l) . These data were obtained for 2D lattice-based 
networks with k = 10“^. We obtain a good collapse of 
the data for / = 0.75 ± 0.05 and u = 2.0 ± 0.1. The 
lower branch, does not vanish as VF ^ oc, due to the 
small, finite value of k. We confirm, however, that for 
W > 1 decreases toward zero as ^ ^ 0 (Supplementary 
Fig. S3). 

In Table [I] we summarize the various values of / and 
f /(/) for different network structures in either 2D or 3D. 
Although the theoretical results above were chosen to 
correspond to connectivity {z) close to the experimen¬ 
tal values for collagen networks, we also studied two very 
different networks—a disordered honeycomb lattice in 2D 
with {z) close to 2 (See Fig. [ 2)3 inset) and a disordered 
FCC lattice with {z) close to (Supplementary Fig. S6). 
Importantly, the near-isostatic case of the FCC lattice at 
{z) 5 exhibits a value f /(f) 1/2, consistent with the 
K ^ scaling reported in Ref. Q for an isostatic net¬ 
work. We note, however, that the individual exponents 
/ and (f> here are only defined for sub-isostatic networks, 
and thus are not ejected to coincide with studies of 
isostatic systems (3, 0|- 

In contrast to most prior work on critical phenom- 


Table I: Critical exponents obtained by simulations of net¬ 
works varying in connectivity, architecture, and dimensional¬ 
ity. 


(^> 

Network 

/ 

f/^ 

2.4 

2d Honeycomb 

0.5 

0.22 

3.2 

2 d phantom 

0.75 

0.36 

3.2 

3d phantom 

0.8 

0.36 

3.6 

2d Mikado 

0.84 

0.4 

5.0 

3d FCC 

1.45 

0.5 


ena for jamming (Hj, rigidity percolation 0,lil and near¬ 
isostatic fibre networks [71, l8||, as well as predictions of 
topological boundary modes in isostatic lattices ( 2 ^ , our 
focus here has been on networks well below the isostatic 
point. This situation is particularly relevant to biol¬ 
ogy, which abounds with structural networks networks 
of biopolymers with connectivity between 3 and 4, well 
below isostaticity in 3D. One challenge in understand¬ 
ing such systems has been their nonlinear mechanical 
response. Recently, a Landau type theory for the non¬ 
linear elasticity of biopolymer gels was proposed using 
an order parameter describing induced nematic order of 
fibres in the gel js^l- Our findings show that the nonlin¬ 
ear mechanics of networks of stiff fibres such as collagen 
can now be understood quantitatively in terms of critical 
phenomena associated with the isostatic point. Impor¬ 
tantly, as we show, there is a line of critical points that 
extends over a wide range of network connectivities, cov¬ 
ering the physiologically relevant range of 3 < z < 4 in 
3D. Moreover, although the linear modulus of collagen 
networks may be finite in this range due to the stabiliz¬ 
ing influence of bending, the nonlinear response can be 
quantitatively captured by the scaling functions in Eqs. 

(HEi). 

Supplementary Information accompanies this paper 
at www.nature.com/naturephysics. 

Acknowledgments We thank M. Vahabi for useful dis¬ 
cussions. This work is part of the research pro¬ 
gramme of the Eoundation for Eundamental Re¬ 
search on Matter (EOM), which is financially sup¬ 
ported by the Netherlands Organisation for Scien¬ 
tific Research (NWO). This work is further sup¬ 
ported by NanoNextNL, a micro and nanotechnol¬ 
ogy programme of the Dutch Government and 130 
partners. 

Author Contributions AS and AJL contributed 
equally to the work. AS, AJL, RR, MS, and 
ECM conceived and developed the model and 
simulations. AS, AJL, RR and MS performed 
the simulations. KAJ and GHK designed the 
experiments. KAJ performed the experiments. All 
authors contributed to the writing of the paper. 













5 


Additional Information The authors declare no 
competing financial interests. Reprints and 
permission information is available online at 
http: / / npg.nature.com/reprintsandpermissions / 

Correspondence should be addressed to GHK or 
FCM. 


Methods Summary 
Network generation 

We model lattice-based networks in 2D and 3D, as well 
as off-lattice (Mikado) networks in 2D. In our lattice- 
based networks, fibres are arranged on a triangular lat¬ 
tice (2D) or a face-centered cubic (FCC) lattice (3D) of 
linear dimension W. In 2D, we randomly select two of 
the three fibres at each vertex on which we form a binary 
cross-link, i.e., enforcing local 4-fold connectivity of the 
network in which the third fibre does not interact with 
the other two. In 3D, where there are 6 fibres crossing 
at a point, we randomly connect three separate pairs of 
fibres at each vertex with binary cross-links to enforce 
local 4-fold connectivity. In both 2D and 3D, the aver¬ 
age connectivity is further reduced below 4 by random 
dilution of bonds with a probability (1 — Pbond)- In order 
to generate 3D networks with connectivity exceeding 4, 
we simply perform random dilution of bonds on a full 
FCC lattice, with an initial connectivity of 12, till the 
desired connectivity is reached. The stretching modulus 
/i and the lattice constant Iq is set to 1 in 2D and 3D. 
Mikado networks are generated by random deposition of 
filaments in a 2D box of size W. A freely hinged cross¬ 
link is inserted at every point of intersection. The depo¬ 
sition continues until the desired average connectivity is 
obtained. 


Collagen Rheology 

For our experiments, rat-tail collagen type I (BD Bio¬ 
sciences, Breda) was polymerized at 37 degrees Cel¬ 
sius, in a in a physiological buffer solution composed of 
DMEMlx solution (diluted from lOx, Sigma) contain¬ 
ing 50 mM HEPES, 1.5 mg/ml sodium bicarbonate, 1% 
EBS (Cibco) and 0.1% antibiotics (pen/strep, Cibco) at 
pH 7.3. Rheology was performed on a stress-controlled 
rheometer (Physica MCR 501; Anton Paar, Craz, Aus¬ 
tria) with a 40 mm cone-plate geometry having an 1 
degree cone angle. A solvent trap was added to main¬ 
tain a humid atmosphere. After 6 hours polymerization, 
stress-stiffening curves were obtained using a differential 
protocol. 
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Concentration scales linearly with the reduced 
bending rigidity k 

The parameter k in our model is naturally related to 
the protein concentration c as follows. For an elastic 
rod of radius r and Young’s modulus fi = 7rr‘^E and 
hz = TTr^EjA, implying that k oc (r//)^ oc since the 
volume fraction (f = 7 rr^^ where p oc 1/P is the total 
fibre length per volume (3|-l5|. Hence, the protein con¬ 
centration c (or p) in experiments can be simply related 
to the reduced bending rigidity k 8iS p k. The theoret¬ 
ical elastic energy involves a summation over all fibres in 
the network and is a function of the strain 7 and the re¬ 
duced bending rigidity k. Moreover, since the modulus K 
involves the energy per unit volume, K is naturally pro¬ 
portional to p. The modulus can therefore be expressed 
as 

K = ppJC ( 7 , k) . (S 6 ) 

In the linear regime, 7 ^ 0, we obtain from the above 
the well known relation K ^ ppk ^ p^. 


Supplementary Information 


Scaling analysis 


Critical strain as the inflection point of 
strain-stiffening curves 


Einite size scaling 


We define the critical strain 7 c as the inflection point 
of the logK vs. log 7 curves, analogous to the determina¬ 
tion of the critical point in a finite size system [l|. Using 
this definition, we can unambiguously extract the criti¬ 
cal strain in a consistent manner for both experiments 
and simulations. In Fig. ISlh. we show the critical strain, 
obtained from simulations on a 2D phantom triangular 
network with {z) 3.2, for different values of the fibre 

bending rigidity. As can be seen, 7 c shows a weak de¬ 
pendence on the k. In the limit of ^ ^ 0, 7 c approaches 
a constant value. Moreover, 7 c also marks the strain at 
which the network exhibits the largest non-afiine fluc¬ 
tuations. Given the displacement field u and the aj/fine 
displacement field of the network, the non-afiine fluc¬ 
tuations can be quantified as (H 


^r(7) 


(II^uNA||2) 


(S5) 


where ST{'y) is referred as differential non-affinity, 
= u — uA is the differential non-afiine displacement 
of a crosslink to an imposed strain dy, I is the typical 
network mesh size and the angular brackets represent a 
network average. In fact, in the limit of ^ > 0, the fluc¬ 
tuations are expected to diverge. In Fig. ISTb . we show 
^F( 7 ) in the neighbourhood of the critical strain. As can 
be seen, the fluctuations grow with decreasing con¬ 
sistent with the idea that the bending rigidity can be 
considered as an auxiliary field. 


The critical behavior in our model can be tested by 
performing finite-size scaling, which is sensitive to the 
divergence of the correlation length. In our system, the 
order parameter K scales as \Aj\^ (Fig. IS^. as system 
size LF ^ oc and A 7 ^ 0. This scaling should be evi¬ 
dent when the correlation length, which scales as 
is much smaller than the system size, where u is the cor¬ 
relation length exponent. At the critical point, when 
the correlation length diverges, the modulus should scale 
with system size as K r\j W such that iF ^ 0 as 
LF ^ oc. This can be summarized in the following scal¬ 
ing relation. 


K = W-f (1A7I. (S 7 ) 

In order to perform finite-size scaling, we chose a finite 
but small bending rigidity. We considered system sizes 
from LF = 40 to a maximum of LF = 200 in 2 D. We chose 
a finite k to avoid the numerical problems associated with 
a rope-like network (^ = 0). In a rope-like network, due 
to finite size effects, the energy density shows a jump at 
the critical strain. Since one needs to take derivatives 
of energy density, it is numerically problematic to un¬ 
ambiguously extract the modulus at the critical strain. 
In Fig. [S^ we show the stiffening curves for k = 10 “^ 
for different system sizes. In Fig. IS^ we present the 
collapse (Eq. dSZl) of stiffening curves obtained for two 
other bending rigidities, k = 10 “^ and 10 “^. As can be 
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Supplementary Figure SI: (a) Critical strain 7 c, obtained from simulations on a 2 D phantom triangular network with 
(z) ~ 3.2, as the inflection point of the logi^ vs. log 7 curves (circles). Squares show the strain at which the non-affine 
fluctuations show a maximum. It is clear that the maximum in the non-affine fluctuations occurs at 7 ~ 7 c. (b) Differential 
non-affinity ^r( 7 ) obtained from the same simulations in (a) as a function of the applied strain. ^r( 7 ) peaks at 7 ~ 7 c. As 
expected, the height of the peak increases with decreasing k since the displacement field of these networks becomes highly 
non-affine as 7 ^ 7 c. 


seen, for the three values of the data collapse is in¬ 
distinguishable except for the lower branch. The lower 
branch converges to a value which scales with k. It fol¬ 
lows that the lower branch extends continuously to zero 
as ^ ^ 0 , consistent with the fact that a central-force 
sub-isostatic network is unstable for 7 < 7 c. 


scaling law. 


|A 7 |^ ^ IA7K V IA7I; 


(S13) 


Equation of state 


It can be seen that for A 7 = 0, the above scaling relation 
correctly reproduces K ^ k^^^ at the critical point. 


In ferromagnetism, the magnetization, m, in pres¬ 
ence of an applied field h and reduced temperature 
t = {T — Tc)/Tc can be captured in the following scaling 
relation Q: 


h / 1//0 \ 

ri ^ / I 1 m \ 


(S 8 ) 


Here, the (±) branch corresponds to t ^ 0 and A and 
P are critical exponents. The analogous quantities for a 
fibre network are the following: 


A 


0 

/3 


/ 

h 

J \ 

k 

w 


A 7 

m 

J \ 

K 

w 


IA 7 K 


(S9) 

(510) 

(511) 

(5 12 ) 


Based on the above analogy, we obtain the following 


In the main text, we use Eq. ()S13D to obtain fit to the 
experimental K ys. 7 data in the following way. We 
first focus on the linear regime. In the linear regime, we 
know from simulations that the modulus (in units of pja) 
scales linearly with k which itself scales as ^ ^ p giving 
rise to the well known (or dependence of the linear 
modulus where c is the protein concentration. It fol¬ 
lows that in order to compare experimental K with that 
obtained from simulations we should first rescale the ex¬ 
perimental K hy so that the rescaled modulus scales 
as ^ c ^ k. As the next step, we obtain the indi¬ 

vidual critical strains, 7 c, for each of the concentrations 
as the inflection point of the \ogK vs. log 7 curve. We 
then consider the experimental data (rescaled by for 
the lowermost concentration along with its 7 c and per¬ 
form fitting to it using Eq. (|S13[) with k as the only free 
parameter. Once we obtain k for a given protein con¬ 
centration, we obtain the k for other concentrations by 
using k ^ c. However, in the fit shown in the main text, 
we have treated ^ as a free parameter. We show in the 
main text that the fit values of k are consistent with the 
expected k r\j c scaling. 
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Supplementary Figure S2: Stiffening curves for 7 > 7 c in a 
limit k ^ 0, K where for this network / ^ 0.75. 

Average connectivity of collagen networks 

In Fig. [S4l we show a typical Scanning Electron Mi¬ 
croscope (SEM) image of a collagen network prepared at 
a temperature T = 37° and concentration of 4 mg/ml. 
Eor SEM, collagen gels (50-100 /il) were polymerized 
overnight inside 5 ml eppendorf tubes in humid condi¬ 
tions. After polymerization, samples were washed three 
times with sodium cacodylate buffer (50 mM cacody- 
late, 150 mM NaCl, pH 7.4) for 30-60 min each, at 
their polymerization temperature. Samples were fixed 
with 2.5% glutaraldehyde in the same buffer for at least 
2 hours. Next, samples were washed three times with 
sodium cacodylate buffer (room temperature) and dehy¬ 
drated with increasing amount of ethanol. After com¬ 
plete dehydration (100 % ethanol), 50 % hexamethyldis- 
ilazane (HMDS) in ethanol was added (under the hood) 
and afterwards replaced after 30 min by 100% HMDS. 
The HMDS was left to evaporate overnight. The samples 
were transported to a stub with carbon tape and sput- 


2D network (phantom triangular, (z) ~ 3.2) shows that in the 


ter coated using a K575X sputter coater (Quorum Tech¬ 
nologies, Gouda, The Netherlands). A layer of 15.4 nm 
of Au/Pd was sputtered using a current of 80 mA. The 
SEM samples were visualized using a Scanning Trans¬ 
mission Electron Microscope (STEM) setup (Verios 460, 
EEI Company, Eindhoven, the Netherlands) using 50 pA, 
5 kV and 4 mm working distance, in immersion mode. 
Two different samples, with number of branches > 100, 
concentration of 4mg/ml and temperature T = 37°, were 
analyzed. The average connectivity of a collagen network 
with concentration 4mg/ml is measured to be 3.3 ± 0.1. 


Critical strain as a function of connectivity 

In Eig. [S5l we show the critical strain as a function 
of average connectivity in a 2D phantom triangular net¬ 
work. In the main text, a schematic version of such a 
figure is shown as Eig. 1. As can be seen in Eig. [SH 7c 
approaches zero as the connectivity approaches the iso- 
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Supplementary Figure S3: (a) Stiffening curves for = 10 ^ for different system sizes, (b) Collapse of the stiffening curves 
according to Eq. (EZI) for different system sizes W. Stiffening curves are obtained for three different bending rigidities (see 
legend). The lower branch converges to a value which scales with the bending rigidity k. In the limit of ^ ^ 0, the lower 
branch will extend all the way down to zero. 


static threshold of (z) = 4 in 2D. In Fig. [SH we mark @| 
as shaded rectangle the connectivity-range relevant for 
collagen networks. In the main text, the overall decrease 
in the critical strain in the experiments (Fig. 3b inset 
in the main text) with increasing concentration of colla¬ 
gen is about 30%. Increasing concentration of collagen is 
expected to result in an increase in the average connec¬ 
tivity of the network. Consistent with this, one can see in 
Fig. [SH that the critical strain decreases with increasing 
connectivity. In fact, jc decreases significantly from 0.2 
to 0.1, as the connectivity increases from 3.2 to 3.6. 

Near isostatic 3D FCC network 

Although the theoretical results shown in the main text 
correspond to connectivity (z) close to the experimental 
values for collagen networks, we also studied a disordered 
FCC lattice with (z) close to 5. The predicted scaling is 
shown in Fig. [SH Such a network is close to the isostatic 
threshold of (z) = 6 in 3D. The near-isostatic case of 
the FCC lattice at (z) ^ 5 exhibits a value f/cp 1/2, 
consistent with the K ^ scaling reported in Ref. Q 
for an isostatic network. The individual exponents are 
/ = 1.45 and (p = 2.9. These critical exponents are 
only defined for sub-isostatic networks, and thus are not 
expected to coincide with studies of isostatic systems (3, 
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Supplementary Figure S4: SEM image of collagen network prepared at a concentration of 4 mg/ml. The branched structure 
of the collagen fibres is visible. The white scale bar represents 2 fim. 
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Supplementary Figure S5: Critical strain as a function of the average connectivity in a 2D network (phantom triangular). 
In the shaded rectangular region, {z) changes from 3.2 to 3.6, a range relevant for collagen. 
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Supplementary Figure S6: Collapse of data from a 3D FCC lattice network with {z) ~ 5. The critical exponents are / = 1.45 
and (j) — 2.9. The ratio f /cj) — 1/2 is consistent with the K ~ scaling reported in Ref. Q for an isostatic network in 3D. 



